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Abstract. 

A computational linear stability analysis of spiral waves in a reaction-diffusion equation is per- 
formed on large disks. As the disk radius R increases, eigenvalue spectra converge to the absolute 
spectrum predicted by Sandstede and Scheel. The convergence rate is consistent with 1/R, ex- 
cept possibly near the edge of the spectrum. Eigenfunctions computed on large disks are compared 
with predicted exponential forms. Away from the edge of the absolute spectrum the agreement is 
excellent, while near the edge computed eigenfunctions deviate from predictions, probably due to 
finite-size effects. In addition to eigenvalues associated with the absolute spectrum, computations 
reveal point eigenvalues. The point eigenvalues and associated eigenfunctions responsible for both 
core and far-field breakup of spiral waves are shown. 
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1. Introduction. Rotating spiral waves are found in many chemical and biolog- 
ical systems and have been the subject of intense study for many years [1011141 [TBI I28j . 
The equations governing these systems are typically of reaction-diffusion type. Al- 
though each system is modeled in detail by specific equations — which are often very 
complex — generic features of the spiral waves can be understood from reaction- 
diffusion equations with simple nonlinearities. Figure 11.11 shows a spiral wave in a 
generic model reaction-diffusion system described in detail in For the model pa- 
rameters in figure n~TI the spiral wave rotates with constant frequency and shape, i.e. 
it is a rotating wave. 




Fig. 1.1. Rotating spiral wave solution of reaction- diffusion equations described in Colors 
indicate the level of the u field with blue used for u near zero and red used for u near 1. The 
wave rotates counterclockwise. The domain radius, R = 80, is approximately 10 times the spiral 
wavelength. Homogeneous Neumann boundary conditions, corresponding to zero chemical flux, are 
imposed at the domain boundary. Model parameters are a = 0.75, b = 0.0006, and e=0.0741 

The focus of our work is a computational study of the linear stability spectra of 
rotating spiral waves such as those shown in figure ITTI To explain the motivation 
behind this study it is necessary to first recall the recent analysis by Sandstede and 
Scheel [211 1221 1281 12~4"| on the spectra of rotating spiral waves. Their work examines 



'Mathematics Institute, University of Warwick, Coventry, CV4 7AL 
twheelerOmaths .Warwick. ac .uk. 
tbarkleySmaths .Warwick. ac .uk. 



2 



P. WHEELER AND D. BARKLEY 



spectra on large bounded disks and on unbounded domains. The results can be 
summarized as follows (see figure [Q)| . On large bounded disks, the linear stability 
spectrum consists of point eigenvalues and what is called the absolute spectrum. The 
absolute spectrum is not actually part of the stability spectrum. However, all but 
possibly a finite number of point eigenvalues converge to the absolute spectrum as 
the domain size tends to infinity. That is, except for finitely many eigenvalues that 
are created through the underlying pattern as a whole, or possibly by the boundary 
conditions, all eigenvalues on large bounded domains are expected to be close to the 
absolute spectrum. The point eigenvalues have well-defined limits as the domain size 
tends to infinity. 

In practice the absolute spectrum must be computed numerically for any given 
reaction-diffusion equation, e.g. |22j . Such computations require discretization in 
only one space dimension and thus are relatively simple compared with computing 
eigenvalues of the full stability problem on a large domain, such as in figure 11.11 

For spiral waves on the unbounded plane, the linear stability spectrum consists 
of point eigenvalues and the essential spectrum. The essential spectrum is continuous 
spectrum and is determined only by the far-field waves trains of the spiral. It too is 
relatively easy to compute numerically in one space dimension. The point eigenvalues 
again depend on the underlying spiral pattern as a whole. 
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Fig. 1.2. Illustration of spectra in the complex plane for spirals on bounded and unbounded 
domains. T, a ts and T, ess represent the absolute and essential spectrum respectively. Points represent 
eigenvalues on a large bounded domain which approach E a j, 3 as the domain size tends to infinity. 
Crosses represent the point spectrum which does not approach £ a f, s as the domain size tends to 
infinity. 

To see how these linear stability spectra may be relevant in practice, we show 
in figure 11.31 simulations of two instabilities of rotating waves on relatively large do- 
mains and the corresponding absolute and essential spectra obtained by Sandstede 
and Scheel [221 ■ ^ n eacn case a rotating wave becomes unstable in a rather dramatic 
fashion and the spiral breaks up. Multiple spiral waves appear in each of these sim- 
ulations shortly after the time shown. In figure ll.3f a) the breakup initiates in the 
central region of the spiral and is referred to as core breakup 1151 122] whereas in 
figure ll.3f b) the breakup first takes place in the outer regions of the spiral and is 
called far-field breakup [i EH1 E21 EH E3 

The case of far-field breakup, figure fOl (right), has been the subject of several 
past studies [Tll2llHl^l ll8H22ll2*5ll26H29| . The breakup can be mostly understood from 
analysis and simulations of one-dimensional systems. While many of these studies are 
based on the complex Ginzburg-Landau equation, results appear to be similar for 
the case of reaction-diffusion equations EI EE]- The typical scenario is that as 
a parameter is varied the spiral first becomes convectively unstable. In a bounded 
domain the onset of convective instability does not generally lead to breakup because 
unstable modes typically propagate away from the core and are not reflected at the 
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Fig. 1.3. Two examples of spiral breakup - core breakup on the left and far-field breakup on the 
right. The top plots show the u chemical field at about the time of breakup in numerical simulations 
in square geometries with homogeneous Neumann boundary conditions. Domains are of length 160 
on a side. The bottom shows the absolute and essential spectrum obtained by Sandstede and Scheel 
for the parameter values used in the simulations. Note that these spectra repeat periodically in 
the imaginary direction, but this can only be seen in (d). Model parameters are: (left) a = 0.75, 
b = 0.0006, and e = 0.0741, (right) a = 0.84, b = -0.045, and e = 0.075. 



boundary. As the parameter is varied further the spiral becomes absolutely unstable. 
Only at the absolute-instability threshold will instability surely occur in a bounded 
domain. Absolute instability corresponds to a growing "global mode" |251 126| . which 
here means an eigenfunction on the bounded domain whose eigenvalue has positive 
real part. 

In the analysis of Sandstede and Scheel, convective instability is signaled by the 
crossing of the essential spectrum into the right half plane |2T][22]- Figure fOT d) shows 
that the spiral is convectively unstable. However, the spiral was already convectively 
unstable prior to the breakup seen in figure 1131 and this is not the cause of breakup. 
The breakup is caused by an eigenvalue with positive real part and the corresponding 
global mode on the finite (bounded) domain. In principle such an eigenvalue could 
be associated with the absolute spectrum or it could be a point eigenvalue. In fig- 
ure 11.3r d) the absolute spectrum is away from the imaginary axis and thus is not 
expected to play a direct role in the far-field breakup. Thus we expect there to be at 
least one positive point eigenvalue not contained in figure rOT d). 

It is worth being clear about potentially confusing terminology. Absolute insta- 
bility is not associated only with the absolute spectrum. The union of the absolute 
spectrum and the point spectrum determine absolute stability. If part of the absolute 
spectrum lies in the right half plane, then the spiral will necessarily be absolutely 
unstable. However, the converse is not true, since absolute instability can arise due 
to point eigenvalues, even if the absolute spectrum lies entirely in the left-half plane. 

We should warn the reader that simulations at the stated parameters in fig- 
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ure ll.3f bl are very sensitive to numerical resolution. As we shall see, this is because 
the particular parameters considered by Sandstede and Scheel are extremely close to 
a transition between far-field and core breakup. 

The case of core breakup, figure rOT leftl. has not been extensively analyzed, in 
large part because one cannot expect to capture much of the spiral core structure in 
one-dimensional studies. (See however Sandstede and Scheel [22] have raised the 
possibility that core breakup may be due to the absolute spectrum. In figure ITTffi c) it 
can be seen that the absolute spectrum is very near the imaginary axis, although it is 
entirely in the left-half plane. (The argument of Sandstede and Scheel is not simply 
that the absolute spectrum is close to the real axis, but the details are not important 
here.) The other possibility for core breakup is that the instability is again due to 
point eigenvalues. The essential spectrum in ll.Sf cl extends into the right-half plane 
and so the spiral is convectively unstable. 

Computing the eigenvalue spectra of spiral waves on large domains has thus be- 
come important. First and foremost, it is desirable to test the absolute spectra of 
Sandstede and Scheel in at least a few cases. The primary issue is whether or not 
eigenvalues tend to the absolute spectrum for typical domains sizes used in studies of 
spiral waves, e.g. domains such as those in figure fOl The theory is still developing 
and we would like to know whether absolute spectra in fact have any implications 
for domains of reasonable size. The other important issue which computations can 
address is the abundance and importance of point eigenvalues not predicted by the 
absolute spectrum. For example, it is desirable to know how many point eigenvalues 
are present within the region shown in ll.3f c1-(d'K how many of these eigenvalues have 
positive real part, and whether or not these are associated with breakup. For these 
reasons we have undertaken the large-scale eigenvalue computation reported here. 

Throughout this paper we shall use point eigenvalue to mean those eigenvalues 
which remain isolated as the domain radius becomes large, as contrasted with the 
eigenvalues associated with the absolute spectrum that approach one another as the 
radius becomes large. To simplify discussion we shall use positive eigenvalue to mean 
an eigenvalue, or a complex-conjugate pair of eigenvalues, with positive real part. 
Similarly we shall use positive eigenfunction to mean an eigenfunction whose corre- 
sponding eigenvalue has positive real part. 

2. Model and Methods. 

2.1. Model. We will consider a standard two-component reaction-diffusion model 5 
given by the equations 



(2.1b) 



(2.1a) 



— = V 2 u + f(u, v) 

dv 2 . . 

— = 6 V v + g(u,v), 



where 



(2.2) 




There is freedom in the choice of g(u, v) and our methods will not depend on this 
choice. However, the results we report will be for the g proposed by Bar and Eiswirth [3] 
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and used by Sandstede and Scheel [221, namely 

-v, 0<u<l/3 
(2.3) g(u,v) = { 1 -6.75u(u- l) 2 -v, 1/3 < u < 1 

L — u, 1 < u 

The equations are posed on a disk of radius R and with homogeneous Neumann 
boundary conditions at r = R: 

^(R,e) = f r (R,9) = 0, 

where r, 9 are standard polar coordinates. For chemically reacting systems these 
are the most natural boundary conditions as they correspond to zero chemical flux 
through the boundary of the domain. Other boundary conditions could give different 
spiral solutions and linear stability spectra on finite domains, but we do not consider 
any other boundary conditions here. 

The parameters of the model are kinetics parameters a, b, and e, and the diffusion 
coefficient 5. If b > the equations model an excitable medium. In this case the 
homogeneous state with u = v = everywhere is linearly stable and finite amplitude 
perturbations are required to initiate waves. The perturbation threshold is set by b/a. 
For b < the equations model an oscillatory medium. In both cases e controls the 
time-scale ratio between the u- and w-equations. We consider e<l corresponding to 
a fast time scale for u relative to v. We shall only report results for the case 6 = 0. 
5 = 1 is the other case commonly considered. As stated at the outset, these equations 
model generic features of spiral waves in a variety of excitable and oscillatory media. 

2.2. Computational tasks. Consider rotating-wave solutions of i|2.1[l rotating 
at frequency u. We use (u* ,v*) to denote such solutions and refer to them as steady 
spirals, since these are steady states when viewed in the frame of reference which 
is rotating with the spiral. Transforming to a system of coordinates co-rotating at 
frequency w, steady spirals obey the equations 

B11* 

(2.4a) o-VV+w— +/(«>*) 

(2.4b) o = ^+ ff ( U >*), 

subject to homogeneous Neumann boundary conditions. These steady-state equations 
can be written in the form 

(2.5) T ( £ ) = 0, 

where T is the nonlinear operator given by the right hand side of l|2.4|l . 

Next, given a steady spiral, we seek the leading part of its linear stability spec- 
trum. Consider the linearized evolution equations, in the rotating frame, for infinites- 
imal perturbations (u',v') of the steady solution (u*,v*): 

Fin* Fin* 
(2.6a) ^ = VV + uj U — + /*(«»' + /„(«»' 

(2.6b) ^ = oj^ + g u (u*,v*)u' + g v (u*,v*)v', 
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where /«,•■■,<?„ denote the derivatives of the kinetic functions. In this frame of 
reference the eigenvalue problem is 

(2.7) C 



where ^ !f J are eigenfunctions, A are the corresponding eigenvalues, and C is 

(2g) c= ( V 2 +uod e + f u (u\v*) f v (u*,v*) 

1 ' ' V g u {u\v*) codg +g v (u*,v*) 

Thus our primary numerical tasks are the solution of steady state problem 1|2.5(1 and 
the determination of the leading eigenvalues of problem Ij2.7|) . 

In addition it is necessary to perform a few simulations of the time-dependent 
equations H2.1(l . e.g., the simulations shown in figure fOl In the case of figure fOl 
the numerical methods are described fully elsewhere I12| and will not be discussed 
here. 

2.3. Computational Methods. Equations l|2.5|l and (|2.7(l are common in large- 
scale numerical bifurcation problems and the computational methods we employ are 
more or less standard 11 . For completeness we provide a basic description of our 
methods and stress a few points concerning implementation which are essential to the 
efficiency of the computations. 

The fields are discretized on a regular polar grid (rj,9k) = (jAr 7 kA8), where 
< j < N r and < k < Ng, plus the center point (0,0). There are thus N r N e + 1 
grid points. The r-derivatives in the differential operators are evaluated using second- 
order finite differences, taking into account the boundary condition at r = R. The 
^-derivatives are evaluated spectrally using Fourier transforms. In this way (|2.5|l and 
(|2.7|) become 

(2.9) F(U*) = 0, 

(2.10) LU = AU 

where the U's are vectors of length N = 2(N r Ng + 1), F is a nonlinear function, and 
L is an N x N matrix. 

Newton's method is used to solve steady state problem l|2.9ll . One iteration of 
Newton's method is 

(2.11) DF(U„)AU„ = -F(U n ) 

(2.12) U n+1 = U n + AU„, 

where DF(U„) is the linearization of F about the current iterate U„. This is the 
same matrix as L except it is evaluated at U n rather than at the steady state U* . 

The work of each Newton's iteration is dominated by the work necessary to solve 
the N x N linear system of equations l|2.11|l for the n th correction AU„. This can 
be done by a direct method if care is taken to order variables to keep the matrix 
bandwidth of DF as small as possible. Let Ujk and Vjk be values at grid point 
(rj,9k). Then these arc ordered in Ui such that the chemical species changes fastest 
with index i, followed by the angular index fc, followed by the radial index j. This 
ordering is not that suggested by (|2.8|l . With the ordering we use the bandwidth is 
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Fig. 2.1. Banded structure of matrix DF or L. Solid black represents full blocks of size 
2Ng X 2Ng mainly due to using spectral representation of 9 -derivatives. The lines are due to the 
second-order finite-difference representation of r- derivatives. 

approximately 4iVg. See figure Even for moderately large discretizations a direct 
method can be used to solve <|2.11[1 . For example on a grid N r x Nq — 600 x 256, N 
is approximately 3 x 10 5 , while the bandwidth is only about 10 3 . 

The only other issue concerning the steady state computations is that the fre- 
quency lu must be found in addition to fields u* and v*. The existence of the addi- 
tional unknown is consistent with the fact that the solution to (|2.4[1 is not unique due 
to the rotational symmetry in 9. One more algebraic equation must be added to l|2.9|> 
to remove the phase degeneracy and thus giving N + 1 equations in N + I unknowns. 
The constraint we add is simply to fix u as some point in the domain. While the 
constraint destroys the banded structure of DF, a Sherman-Morrison technique ^5] 
is used to find solutions of the augmented linear system using only the banded DF. 

We now describe our computations of the leading eigenvalue spectrum of L. The 
basis of our approach is to employ a Cayley transformation to transform the eigenval- 
ues we seek (those with largest real part) to dominant eigenvalues (of largest magni- 
tude), and then to find iteratively dominant eigenvalues of the transformed operator. 
Reference )17| gives a nice review of such methods. Specifically, we consider the matrix 
A defined by the Cayley transform 



where £ and rj are real parameters and I is the identity. Letting the ji and A be the 
eigenvalues of A and L respectively, we have the relation 



The parameters £ and rj can be adjusted so as to map the regions of interest in the 
A-plane to large magnitude in the //-plane. Using the predicted absolute spectra of 
Sandstede and Scheel it is easy to find appropriate values of £ and rj. Figure 12.21 
shows the effect of the Cayley transform on the absolute spectrum for one of the cases 
predicted by Sandstede and Scheel [22 (the other case is similar) for the values of £ 
and rj used in our computations: £ = —0.4 and rj — 4.0. In the /i-plane we include 
transforms of two periodic repeats of the absolute spectrum, one corresponding to 
the larger imaginary part and one to smaller imaginary part (in the A-planc ) . These 
repeats are outside the region of the A-plane shown. Most of the eigenvalues of L lie 
far to the left in the A-plane (outside the range of the figure). These are all mapped 
to near the origin by l|2.13[l . 

There is no need to form or store the matrix A in order to iteratively calculate 
its eigenvalues. All one needs is the ability to compute AU for arbitrary U. This 



(2.13) 



A=(a + L)- 1 (r ? I + L), 
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Fig. 2.2. Plot showing the effect of the Cayley transform on the absolute spectrum in \l. !W c ) . 
Left is the original and right is transformed spectrum. The right includes two periodic repeats of the 
absolute spectrum which are outside of the region shown on the left. Shading indicates magnitude 
of eigenvalues after the transformed. 



is accomplished using the same basic technique as in Newton's method. Letting 
U' = AU we see that U' obeys 

(2.15) (£I + L)U' = (»?I + L)U. 

However, (£1 + L) has the same structure, in particularly the same bandwidth, as L 
and (77I + L) requires mostly the same computations as evaluating F. Thus we act 
with A on U by computing (771 + L)U to form a right-hand side and then solving a 
linear system with matrix (£1 + L). Since this is a fixed matrix, for any given L, it is 
factored only once for all actions of A. 

Dominant eigenvalues A are easily found by subspace iteration |171 I27| . This 
method guarantees that we can obtain any required number of dominant eigenvalues 
to arbitrarily high precision. While Arnoldi's method generally converges faster, in 
practice we find that with this method all eigenvalues we require are not found with 
high enough precision. While there are methods, such as block Arnoldi, which could 
probably address this, we have used subspace iteration. From the eigenvalues fi of A 
we invert (|2.14|) to find the required eigenvalues A. 

2.4. Accuracy. We conclude this section by considering the accuracy of our 
computations and providing details of numerical parameters used in the results re- 
ported. The sources of error are the following: 

1. Discretization error of the steady state problem, i.e. approximation of (|2.5[1 
by (22). 

2. Residual error arising from determining the roots of (|2.9(l to a finite accuracy. 

3. Discretization error of the eigenvalue problem, i.e. approximation of (|2.7|) by 
(I2TTH) . 

4. Residual error arising from computing eigenvaluc/cigcnfunction pairs l|2.10|l 
to finite accuracy. 

The two residual errors are least important. We always iterate sufficiently to 
reduce these to negligible size. The following hold for all reported results. Solutions 
U* of (22) are found such that ||F(U*)|| < 1(T 8 . Solutions of (2H5) are found such 
that || LU - AU|| < 1CT 8 , where ||U|| = 1. The norm is the vector 2-norm. The 
dimension k of the subspaces used in subspace iteration are: k — 30 for R = 20 and 
R = 40, and k = 75 for R = 80. In all cases we stop iterations when ~ 0.7fc of the 
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Fig. 2.3. Graphs illustrating the convergence of steady states (top) and eigenvalues (bottom) 
as function of grid resolution Ar 2 , where ui is the spiral freguency and A c is a complex eigenvalue. 
The domain radius is R = 40. Crosses are with Ng = 256. For Ar = 0.1333 computations have been 
repeated with Ng = 128 are results are shown with circles. Parameters are a = 0.75, b = 0.0006, 
and e = 0.0741. 

eigenvalue-eigenvector pairs have residual less than 10~ 8 . In the case of R = 80, we 
thus obtain 53 pairs with the required residual. We initially start with a subspace 
generated from k random vectors, but we restart iterations from previous runs when 
necessary. 

The discretization errors are dominated, in both the steady state and eigenvalue 
computations, by the second-order finite-derivative approximation to the r-derivatives 
in the Laplacian operator. This is expected since the ^-derivatives are computed with 
spectral accuracy. Thus the dominant error in the computations depends on Ar in 
a well-understood way. Figure 12.31 shows examples of how results from steady state 
and eigenvalue computations scale with Ar 2 . The domain radius is R = 40, half the 
maximum considered in our work. For the steady states we show the frequency ui 
and for the eigenvalues we show the magnitude of A c , a complex eigenvalue associated 
with core breakup that will be considered in detail later in the paper. In both cases 
the second-order convergence is evident. We are interested only in leading eigenvalues 
(roughly 10 2 out of 10 5 ) all of which correspond to eigenfunctions with variation on 
approximately the same spatial scale (roughly the wavelength of the underlying spiral) 
so that the effects of the finite-difference discretization will be approximately the same 
for all eigenvalues we report. 

Based on these plots, we use Ar = 0.1333 for all results reported in § At 
this value of Ar we have performed computations at both Ng = 128 and Ng = 256. 
These results show clearly that Ng = 128 gives sufficient resolution for domain radius 
R = 40. Therefore, for radii up to at least R = 80, Ng = 256 should produce smaller 
errors than the already small errors due to the radial discretization. In summary, for 
all results in §[3]we use Ar = 0.1333 and Ng = 256. 

3. Results. 

3.1. Spectra. We begin with results for the eigenvalue spectra. Figure lTTI shows 
leading eigenvalues of L for the two cases considered by Sandstede and Schcel. The 
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spectrum corresponding to core breakup is at the top and the spectrum corresponding 
to far-field breakup is at the bottom. In each figure eigenvalues computed for three 
domain radii, R = 20, R = 40, and R = 80, are plotted as points with dashed lines 
connecting eigenvalues associated with the absolute spectrum. For comparison, the 
absolute and essential spectra obtained by Sandstede and Scheel are shown as solid 
curves. 
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Fig. 3.1. Eigenvalue spectra, (a) Spectrum corresponding to core breakup and (b) spectrum 
corresponding to far-field breakup. In each case eigenvalues are shown for three domain radii: 
R = 20 (green diamonds), R = 40 (red crosses), and R = 80 (blue circles). For each radius, the 
eigenvalues associated with the absolute spectrum are connected with lines. Point eigenvalue are not. 
Parameters for (a): a=0.75, b=0.0006, e=0.0741, u=1.71. Parameters for (b): a=0.84, b=-0.045, 
e=0.0751, u=1.50. 
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Before considering either case in detail, we note that the predominant feature 
in both is the presence of many eigenvalues which approach the predicted absolute 
spectra as the domain radius increases. Each doubling of R results in approximately 
a halving of the distance of eigenvalues to the absolute spectrum, thus supporting 
l/R convergence to the absolute spectrum. Furthermore the density of eigenvalues 
approximately doubles with each doubling of R. We return to this while considering 
each case in more detail. It should be noted that in the far-field case we have not 
obtained all eigenvalues associated with the periodic repeats of the absolute spectrum 
for R — 80 due to the difficulties of computing these with sufficient accuracy. This 
results in an abrupt termination of the eigenvalue branches for R — 80 at the top and 
bottom of figure mT bl. 

Consider first the spectrum in figure l3~l7 a) corresponding to core breakup. Within 
the region of the complex plane shown, there are three point eigenvalues. All other 
eigenvalues are associated with the absolute spectrum. Specifically, we find three 
eigenvalues which are insensitive to the domain radius and which are separated from 
the absolute spectrum. Of these, one is the zero eigenvalue arising due to rotational 
symmetry. There are three points indistinguishable from zero in figure UTTT a) corre- 
sponding to the three domain radii studied. The other two point eigenvalues are a 
complex-conjugate pair at approximately 0.050 ± 0.543i. As we shall show in § 13.31 
these eigenvalues are associated with core breakup. We will denote them by A c . 
(These are the eigenvalues considered in the convergence study in figure l2~3"l ^ All 
other eigenvalues vary with domain radius and approach the absolute spectrum as 
the radius becomes large. 

The enlargement in figure mT a - ) clarifies the situation around the complex point 
eigenvalues. Even on the enlarged scale the point eigenvalues for R = 40 and R = 80 
coincide. At R = 20, however, there are two nearby eigenvalues. The lower one cor- 
responds to the absolute spectrum (indicated by the connecting dashed line) because 
it moves, as R is increased, toward the absolute spectrum. The other eigenvalue con- 
verges, as the R is increased, to the point eigenvalue. Note that while the edge of the 
eigenvalue branch associated with the absolute spectrum is near the point eigenvalue 
at R = 20, the branch does not approach the point eigenvalue as the domain becomes 
small. It is nevertheless interesting that the point eigenvalue is near the edge of the 
absolute spectrum. We find this throughout and return to this in the conclusion. 

As already noted, the distance of eigenvalues to the absolute spectrum is approx- 
imately proportional to l/R and the density of eigenvalues is approximately propor- 
tional to R. Because we are not able to extend the computations significantly beyond 
the radius R = 80 (already quite large) there is not enough data to draw strong con- 
clusions about the these scalings. In particular it is not clear from the data whether 
or not the scaling in the vicinity of the edge of the absolute spectrum is different from 
elsewhere. Near the edge of the spectrum the eigenvalues are more dense and closer 
to the absolute spectrum than elsewhere. More importantly we observe a curving 
and perhaps folding, at the edge of the eigenvalue branch as the radius becomes large. 
This would again suggest a different scaling at the spectrum's edge, but the numerical 
results are inconclusive. 

We have focused our study on the eigenvalues within the region shown in fig- 
ure E^^a), but we have computed some eigenvalues out side of this region. In partic- 
ular our iterative technique frequently finds eigenvalues associated with the periodic 
repetition of the absolute spectrum in the complex plane. We have not attempted to 
resolve the details of the other eigenvalue branches. Also there is a complex-conjugate 
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pair of point eigenvalues near ± iuj due to approximate translational symmetry. 

Now consider the spectrum corresponding to far-field breakup. In figure 13. If f)) 
we find four complex-conjugate eigenvalue pairs that we can clearly classify as point 
eigenvalues. One of these pairs has small positive real part and hence the spiral 
wave is absolutely linearly unstable, see section § 13.31 Again we observe that the 
point eigenvalues, except for the zero eigenvalue, appear near the edge of the absolute 
spectrum. 

We observe approximately 1/R convergence of eigenvalues to the absolute spec- 
trum. The only apparent deviation is again at the edge of the spectrum. In this case 
we do not observe any curving of the eigenvalue branch seen in the enlargement in 
figure mT a): however, we do find that the right-most point of the computed branch 
does not appear to approach the edge of the predicted absolute spectrum. One pos- 
sibility is that this last eigenvalue is in fact a point eigenvalue very close to the edge 
of the absolute spectrum. 

3.2. Eigenfunctions. In figures l3~2"l and l3~3l we plot eigenfunctions for a repre- 
sentative selection of eigenvalues. All eigenfunctions have been obtained on a domain 
with R — 80, the largest we consider. Each eigenfunction is shown in two formats. In 
the left column eigenfunctions are visualized on the computational domain. Specifi- 
cally the -u-field of the real part of the eigenfunction is plotted with black used where 
the field values are near zero. The imaginary parts of the eigenfunctions and not fun- 
damentally different. In the right column the modulus of each eigenfunction is shown 
as a function of radius. Specifically, we generate 16 radial sections \u(r, 9 S )\ where 
9 S = stt/8 and s = 0, 1, . . . , 15 and plot all 16 sections simultaneously as a function of 
r. The envelope of these sections gives a simple representation of the modulus of the 
eigenfunction as a function of r. 

The technique used by Sandstede and Scheel [231 to obtain absolute and essen- 
tial spectra also predict large-r behavior of eigenfunctions. The main prediction is 
that if an eigenvalue is to the left of the essential spectrum then the corresponding 
eigenfunctions will be exponentially growing at large r, whereas if an eigenvalue is 
to the right of the essential spectrum then the corresponding eigenfunctions will be 
exponentially decaying at large r. In addition to the general prediction, the numerical 
technique employed by Sandstede and Scheel to obtain spectra for specific problems 
also provides the growth/decay rates for eigenfunctions. These rates [20] are indicated 
by the (red) lines in the right column of figures 13.21 and 13.31 For eigenfunctions corre- 
sponding to the absolute spectrum, the predicted exponential growth rates have been 
taken from the point on the absolute spectrum nearest to the computed eigenvalue. 
Only the slope of the lines is relevant. The intercept is chosen for ease of comparison 
with the eigenfunctions. 

Consider first the eigenfunctions in figure 13.21 corresponding to the case of core 
breakup. The top eigenfunction is the zero mode due to rotational symmetry. This 
eigenfunction is given by the ^-derivative of the underlying spiral wave and hence 
closely resembles the spiral. The eigenfunction neither grows nor decays at large r. 

Figure GOf b) shows the eigenfunction corresponding to the positive complex- 
conjugate point eigenvalues A c . Since the eigenvalues are to the right of the essential 
spectrum the eigenfunction decays with r. While there is generally good agreement 
between the computed eigenfunction and prediction, there is some deviation from 
prediction that is more pronounced at larger r. 

Figures l3.2f cV(e') show three eigenfunctions associated with the absolute spec- 
trum. The agreement between the computed eigenfunctions and predictions is ex- 
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Fig. 3.2. Representative eigenf unctions for parameters corresponding to core breakup. The left 
column shows the real part of the each eigenf unction. The u-field is plotted with black used where 
the field is near zero. Right column shows the r dependence of \u\ with predicted growth/decay rate 
also shown with lines (see text). Plot at the far right is a guide to the corresponding eigenvalues, 
(a) zero (rotational) eigenvalue, (b) positive eigenvalue X c , and (c), (d), (e) three representative 
eigenvalues associated with the absolute spectrum. Parameters as in fiaure T3~l\f a). 

tremely good away from the edge of the absolute spectrum. Near the edge the agree- 
ment is less good. In particular, eigcnfunctions arc not pure exponential, even at large 
r, and the computed eigcnfunctions systematically grow more slowly than prediction. 
While not shown, we find that the eigenfunctions computed on smaller domains and 
show even slower growth as a function of r. This would suggest that the deviation 
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Fig. 3.3. Representative eigenf unctions for parameters corresponding to far-field breakup. Same 
format is used as in figure \3.'A The rotational eigenfunction is not shown. Two eigenf unctions 
corresponding to point eigenvalues are shown with (a) slightly positive eigenvalue and (b) slightly 
negative eigenvalue, (c), (d), (e) are three representative eigenvalues associated with the absolute 
spectrum. Parameters as in fiaure \3.lV b). 

shown in figure l3~^T c) is due to finite domain size. 

Figure 13.31 shows eigenfunctions corresponding to parameters for which far-field 
breakup occurs. We plot eigenfunctions corresponding to two of the complex-conjugate 
point eigenvalues and show three representative eigenfunction associated with the ab- 
solute spectrum. The eigenfunction corresponding to the zero eigenvalue is similar to 
figure EOT a) and is not shown. 
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The eigenfunctions corresponding to the point eigenvalues agree very well with 
prediction. The growth rate of the positive eigenfunction in figure I3~37 a) is very small 
since the corresponding eigenvalue is almost exactly on the essential spectrum (fig- 
ure GO^b)). This is a fortuitous situation which illustrates nicely that the essential 
spectrum delimits the crossover from growth to decay of eigenfunctions. While quan- 
titatively the agreement is very good, there is a qualitative disagreement between the 
computed eigenfunction and prediction. Our computed eigenfunction is growing with 
r, indicating that the eigenvalue is actually slightly to the left of the essential spec- 
trum, whereas in figure l3~TT b) the eigenvalue it slightly to the right of the essential 
spectrum and the predicted exponent is slightly negative. The quantitatively the dif- 
ference is very small and is likely due to a small numerical difference, e.g. a difference 
in the value of to found in our computations and that used by Sandstede and Scheel. 
The closeness of this eigenvalue to the essential spectrum is just by chance. If parame- 
ters are changed the eigenvalue moves away from the essential spectrum. It is because 
of this closeness to the eigenvalue to essential spectrum that numerical simulations at 
these parameter values are so sensitive to numerical resolution (as noted in §[T|l. 

The eigenfunction in figure l3~37 b^ is exponentially growing since the corresponding 
eigenvalue is to the left of the essential spectrum. There are no observable deviations 
from pure exponential growth at large r and the agreement with prediction is very 
good. 

The three eigenfunctions associated with the absolute spectrum show the same 
trend as in figure 13.21 The agreement between the computed eigenfunctions and 
predictions is extremely good away from the edge of the absolute spectrum while 
near the edge eigenfunctions are not pure exponential and systematically grow more 
slowly than prediction. This case is even more striking than figurc l3~2l for the following 
reasons. The growth rates in figure l3~2l are much larger than figure l3~31 (Note the scale 
change.) The numerical values representing the eigenfunctions span a larger range and 
yet the computed eigenfunctions away from the edge agree very well with predictions. 
There is every reason to believe that these eigenfunctions are numerically well resolved 
within the finite domain. Unlike the case in figure l3~2T c). here the eigenfunction closest 
to the edge of the spectrum, figure l3~ffi c). deviates from exponential growth only at 
large r. There is a clear range r, up to approximately r = 40, where the eigenfunction 
agrees well the predicted exponential growth. This strongly suggests that the lack 
of agreement is due to finite-size effects. It is nevertheless interesting that these are 
more pronounced near the edge of the spectrum. 

3.3. Implications for breakup. We now return to the issue of spiral breakup 
discussed at the outset (figure ll.3|l . We begin with the case of core breakup. Recall 
that while this was treated by Sandstede and Scheel |22l |2~3"| , they were not able to 
draw definite conclusions about the role of the absolute and essential spectrum in core 
breakup. It is already clear from the spectra in figure 13.11 that the steady spiral is 
linearly unstable due to the presence of positive point eigenvalues A c . Here we present 
additional nonlinear simulations of the breakup. 

Figurc rOl shows the time evolution from two different initial conditions each com- 
posed of the steady spiral plus a small amount of one of the computed eigenfunctions. 
The amplitude plotted is defined as A = ming \\U* — ReU\\, where Re is a rotation by 
angle 9. The norm is the vector 2-norm. Essentially A is the norm of the difference 
between the u-field of the steady spiral, U*, and the it-field of the nonlinear solution 
U. The minimization over rotation is included to take into account the small drift of 
the nonlinear solution relative to the steady spiral. 
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Fig. 3.4. Two time series showing the evolution starting from different perturbations of the 
steady spiral with parameters leading to core breakup. In one case (solid green) the perturbation is 
the eigenfunction corresponding to positive eigenvalue A c [(b) in fiaure iy.SH . In the other (dashed 
red) the perturbation is the eigenfunction corresponding to right-most eigenvalue associated to the 
absolute spectrum [(c) in Haure VJ ."M . The dotted blue curve shows exponential growth at rate given 
by A c . Parameters are as in fiaure \S.lV a). 

Consider the evolution starting from the initial condition formed from the posi- 
tive eigenfunction corresponding to A c . Accompanying visualizations are presented in 
figure 1X51 The dynamics is initially linear, obeying the exponential growth dictated 
by the real part of A c . After a short time the growth becomes nonlinear and almost 
immediately core breakup occurs [figure rT^T c): time 25]. Beyond this time the am- 
plitude A loses most of its meaning. Visualizations at much later times are shown. 
One of the more striking aspects of the breakup is that it occurs at r ~ 20, not at the 
center of the spiral. This radius is near where the unstable eigenfunction has maximal 
magnitude. Visually one sees the similarity between the nonlinear breakup and the 
unstable eigenfunction in figure rTKl 

The initial nonlinear growth in figure[!0]is faster than linear. This means that, at 
lowest order, the effect of nonlinearity on the instability is destabilizing. Such behavior 
occurs, for example, sufficiently close to a subcritical bifurcation, e.g. ^Hj- This 
nonlinear destabilization is consistent with the fact that small positive eigenvalues 
lead to the dramatic breakup of the spiral wave. If nonlinearity were stabilizing, one 
would expect the linear instability to saturate in a state resembling a superposition 
of the original spiral and a small amount of the unstable eigenmode (as occurs for 
example in spiral meandering OH IB])- We note that in systems such as this one the 
behavior can change very rapidly with parameters following a bifurcation |25U26j . and 
hence we are not able to conclude that the nonlinear growth follows from a subcritical 
bifurcation, only that at these parameter values it is destabilizing. 

For completeness we have also computed the nonlinear evolution from an initial 
condition formed from the eigenfunction corresponding to the right-most eigenvalue 
of the absolute spectrum [(c) in figure EO] . Figure |^| shows the initial dynamics from 
this simulation. Not surprisingly A does not change much on the scale of figure 
because the associated eigenvalue is very near zero. The simulation eventually leads 
to core breakup if run long enough. However, this is simply because the steady 
spiral is linearly unstable. When breakup does eventually occur, it follows the same 
route (e.g. same exponential growth) as for the initial condition based on the positive 
eigenfunction. 

The conclusion is that, in this case, core breakup is due to nonlinear effects 
following from linear instability due to a complex-conjugate pair of point eigenvalues. 
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Fig. 3.5. Snapshots of evolution from the perturbed steady spiral leading to far-field breakup 
(solid green curve in figure The u-field is shown with u ~ blue and u ~ 1 red. (a) t = 0, 

(b) t = 20 (~ 5 periods), Jcjt = 25, fcij t = 85 (~ 23 periods), (e) t=120 (~ 33 periods), (}) 
eigenfunction. Parameters are as in figure IXTK a ) . 



The absolute spectrum plays no direct role in the spiral breakup. 

Next we briefly consider far-field breakup. We have directly computed the eigen- 
function associated with absolute instability causing far- field breakup, figure ETffi a). 
The leading eigenfunction shows exactly the long wavelength modulation of the steady 
spiral expected for this instability P |1 M CHI \M Figures and 1X71 show the 
dynamics following from the steady spiral perturbed by the unstable eigenfunction. 
The dynamics are initially that of exponential growth with the expected growth rate. 
The growth becomes nonlinear and long-wavelength modulation of spiral becomes 
visible (time 80 in figure ETTj) . Shortly thereafter the spiral breaks near the domain 
boundary. At these parameter values the eigenfunction grows weakly with radius, as 
seen in figure EOT ah and for this reason one would expect a preference for breakup 
near the domain boundary. However, in this case the qualitative character of eigen- 
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Fig. 3.6. Two time series showing the evolution starting from different perturbations of the 
steady spiral with parameters leading to far-field breakup. In one case (solid green) the perturbation 
is the eigenfunction corresponding to positive eigenvalue [(a) in figure 13., yt . In the other (dashed 
red) the perturbation is the eigenfunction corresponding to the weakly stable point eigenvalue [(b) in 
figure [Xffi . The dotted blue curve shows exponential growth at rate given by the positive eigenvalue. 
Parameters are as in fjgureUlYb). 

function depends sensitively on parameters and for slightly different parameters the 
eigenfunction may decay weakly with radius. 

The far-field case is similar to the core breakup case in most other respects. The 
nonlinear growth in figure 13.61 is faster than linear. No other eigenvalues appear 
to play a direct role in the far-field breakup. Figure 13.61 shows the evolution from 
an initial condition formed from the eigenfunction corresponding to the complex- 
conjugate point eigenvalues near the imaginary axis [(b) in figure HO] . 

4. Conclusions. In this paper we have examined in detail the linear stability 
spectra and associated eigenfunctions for spiral waves in large domains. Everywhere, 
except possibly near the edges of the absolute spectra, numerically computed eigenval- 
ues and eigenfunctions agree extremely well with the results of Sandstede and Scheel. 
Our results answer the question posed at the outset. Absolute spectra can be rele- 
vant and predictive for typical domain sizes used in studies of spiral waves. Even in 
domains containing only a few spiral wavelengths (the case R = 20) eigenvalues show 
signs of the absolute spectrum - they lie along curves located roughly in the correct 
part of the complex plane. For domains containing five spiral wavelengths or more 
(R > 40) eigenvalues lie quite close to the absolute spectra. Of course these results 
are for the particular equations and parameters studied here and absolute spectra 
will not necessarily be such good predictors for domains of these sizes in other sys- 
tems. Nevertheless, in at least two cases, one with excitable dynamics and one with 
oscillatory dynamics, absolute spectra are predictive. 

The computed eigenvalues support convergence to the absolute spectrum inversely 
with the domain radius, at least away from the edge of the absolute spectrum. In 
most cases eigenfunctions agree with the exponential forms deduced by Sandstede 
and Scheel. This is even the case for point eigenvalues not associated with the abso- 
lute spectrum. Near the edges of the absolute spectrum, however, eigenfunctions do 
not exhibit exponential growth at large radius, even in the largest domains we have 
considered. While results suggest that this is due to finite-size effects, more work is 
necessary to understand the behavior of eigenvalues and eigenfunction near the edges 
of the spectrum. 

We have computed the positive point eigenvalues giving rise to both core and 
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Fig. 3.7. Snapshots of evolution from the perturbed steady spiral leading to far-field breakup 
(solid green curve in figure The u-field is shown, (a) t = 0, (b) t = 50 f~ 12 periods), (c) 

t = 80 f~ 19 periods), (d) t = 85, (e) t=120 (~ 29 periods), (f) eigenf unction. Parameters are as 
in fiaure WlY b). 



far-field breakup of spiral waves. The essential difference between the two cases is the 
form of the eigenfunctions. For core breakup the eigenfunction has a maximum not 
far from the core of the spiral and decays at large radius. For far-field breakup the 
eigenfunction grows with radius. Moreover, the far- field eigenfunction shows long- 
wavelength modulation known to precede far-field breakup. Nonlinearity also plays a 
role in breakup and we have presented nonlinear simulations showing the destabilizing 
effects of nonlinearity in each case. 

The most intriguing aspect of this work is that all point eigenvalues we have found 
appear near edges of the absolute spectra. This may be a coincidence, but it would 
not seem so from figure EH1 We leave this as an open problem. 
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